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We present improvements of a recently introduced numerical method [Arrigoni et al, Phys. Rev. 
Lett. 110 , 086403 (2013)] to compute steady state properties of strongly correlated electronic 
systems out of equilibrium. The method can be considered as a non-equilibrium generalization of 
exact diagonalization based dynamical mean-held theory (DMFT). The key modihcatlon for the non¬ 
equilibrium situation consists in addressing the DMFT impurity problem within an auxiliary system 
consisting of the correlated Impurity, Ns uncorrelated bath sites and two Markovian environments 
(sink and reservoir). Algorithmic improvements in the impurity solver allow to treat efficiently 
larger values of Ns than previously in DMFT. This increases the accuracy of the results and is 
crucial for a correct description of the physical behavior of the system in the relevant parameter 
range including a semi-quantitative description of the Kondo regime. To illustrate the approach we 
consider a monoatomic layer of correlated orbitals, described by the single-band Hubbard model, 
attached to two metallic leads. The non-equilibrium situation is driven by a bias-voltage applied to 
the leads. For this system, we investigate the spectral function and the steady state current-voltage 
characteristics in the weakly as well as in the strongly interacting limit. In particular we Investigate 
the non-equilibrium behavior of quasi-particle excitations within the Mott gap of the correlated 
layer. We find for low bias voltage Kondo like behavior in the vicinity of the insulating phase. In 
particular we observe a splitting of the Kondo resonance as a function of the bias voltage. 

PACS numbers: 71.27.-|-a 47.70.Nd 73.40.-c 05.60.Gg 


I. INTRODUCTION 

The recent impressive experimental progress in tailor¬ 
ing different microscopically controlled quantum objects 
has prompted increasing interest in correlated systems 
out of equilibrium. Of particular importance are corre- 
late d h eterostructurej^^, quantum wire^ and quantum 
dot^^ with atomic resolution, experiments in ultra cold 
atomic gases in optical lattices^^ff^, as well as ultrafast 
laser spectroscopjff^ff^. 

The theoretical description and understanding of these 
experiments in particular and of complex strongly cor¬ 
related systems in general presents major challenges to 
theoretical solid state physics. For this purpose differ¬ 
ent theoretical approaches have been developed. For the 
equilibrium situation one of the most powerful methods 
is dynamical mean-field theory (DMFTjfisHlII^ which is a 
comprehensive, thermodynamically consistent and non- 
perturbative scheme. The only approximation in DMFT 
is the locality of the self-energy, which becomes exact in 
infinite dimensions, but usually it is a good approxima¬ 
tion for two and three spatial dimensions. The key point 
of DMFT is to map the original problem onto a single 
impurity Anderson model (SIAM)I^ whose parameters 
are determined self-consistently. For this purpose, sev¬ 
eral classes of so-called impurity solvers were developed. 
Among them, the most powerful methods are the numer¬ 
ical renormalization grou p (NR G) approadPsHUl^ Quan¬ 
tum Monte Carlo (QMC)pSlS2l and exact diagonalization 

(ED^pm. 

Prompted by the success of DMFT for equilibrium sys¬ 
tems, the approach was extendecP^HUl to deal with time- 
dependent problems within the nonequilibrium Green’s 



FIG. 1. (Color online) Schematic representation of the sys¬ 
tem, consisting of the correlated layer (red) with local Hub¬ 
bard interaction U and on-site energy £c, sandwiched between 
two semi-infinite metallic leads (blue), with on-site energies ei 
and Er respectively. The hopping between neighboring sites 
of the correlated layer is tc, while the one for the left (right) 
lead is ti (G). Hybridization between the left (right) lead and 
the correlated layer is vi (vr)- A bias voltage ^ — pi — fir is 
applied between the leads. 


function approach originating from the works of KubcP, 
Schwinger®, Kadanoff, BayirP® ^nd KeldysiP. Simi¬ 
lar to the equilibrium case, also non-equilibrium DMFT is 
based on the solution of an appropriate (non-equilibrium) 
SIAM. Despite the fact that many approaches have 
been suggested to solve such impurity problems (see, 
e.g. Refs. [351 [33 |35J |3SJ not all of them are 

suited for non-equilibrium DMFT. In addition, many of 
these are only reliable for short times and cannot treat 
long time behavior and accurately describe the steady 
state. Therefore, developing a non-perturbative impurity 
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solver, which can treat reliably the steady state behavior 
of the SIAM is quite a challenge. 

The non-equilibrium approach, that will be presented 
in this paper, has its root in the exact diagonalization 
(ED)-based DMFT (ED-DMFT). In equilibrium ED- 
DMFT one replaces the infinite bath by an auxiliary fi¬ 
nite non-interacting electronic chain whose parameters 
are determined by a fit to the DMFT hybridization func¬ 
tion A. This cannot be trivially extended to the steady 
state situation. First of all, due to the fact that the auxil¬ 
iary system is finite, there is no dissipation and a proper 
steady state is never reached. An additional technical 
aspect is that the spectrum of the auxiliary system is 
discrete and therefore, the fit in real frequencies is prob¬ 
lematic. But only in the equilibrium case one can cir¬ 
cumvent this problem by introducing a fit in Matsubara 
space .1^ A possible solution to these problems was sug¬ 
gested by us in Refs. ICT] and [5^ with an approach which 
enables direct access to the steady state properties of 
the correlated impurity problem. The basic idea is that 
in addition to a finite number of bath sites coupled to 
the impurity, as in equilibrium ED, two Markovian envi¬ 
ronments are introduced, which act as particle sink and 
reservoir. This auxiliary model represents an open quan¬ 
tum system with dissipative dynamics, which allows to 
properly describe steady state situations. The behav¬ 
ior of this auxiliary non-equilibrium impurity problem 
is described by a Lindblad master equation, which can 
be solved exactly by numerical approaches such as full 
diagonalizatioiPil, non-Hermitian Krylov spac^^ or ma¬ 
trix product state (MPS) methods.^ Its solution allows 
to determine both, the retarded and the Keldysh self¬ 
energies, which are required by the DMFT loop, with 
high accuracy. Here, in particular, we apply the Krylov 
space approach of Ref. |62] to solve the DMFT impu¬ 
rity problem. This yields a much better accuracy than 
in Ref. which allows us to resolve the splitting of the 
quasi-particle resonance as a function of the bias voltage. 

The paper is organized as follows: In Sec. |II A| we 
shortly introduce the Hamiltonian of the system, while 
in Secs. |HB| and |HC| we give an overview over steady- 
state DMFT within the non-equilibrium Green’s function 
formalism. In Sec. |HD| we discuss the auxiliary master 
equation approach, with focus on details of our imple¬ 
mentation. Afterwards in Sec. |HI| we present our results 
for a simple correlated interface. In particular, in Sec. 
|HI A| we benchmark the accuracy, while in Secs. mE 
and HIC| the steady state current and spectral functions 
are investigated, respectively. Finally in Sec. |IV| we give 
concluding remarks and an outlook. 


II. MODEL AND METHOD 


Fig .[^, which consists of a correlated infinite and transi¬ 
tionally invariant layer (c), with local Hubbard interac¬ 
tion U, on-site energy = —U/2 and nearest-neighbor 
hopping amplitude tc, sandwiched between two semi¬ 
infinite metallic leads (a = l,r), with on-site energies Eq, 
and nearest-neighbor hopping amplitudes ta. The leads 
are semi-infinite and translationally invariant in the xy 
plane (parallel to the correlated layer). The hybridiza¬ 
tion between lead a and the correlated layer is Va (See 
Fig. 0. A bias voltage d) is applied between the leads. 
The Hamiltonian reads 


H = He + ^ + Hcoup . (1) 

OL—l^r 

Here 

Tic = -tc ^ 4cr<^ja + ^ X! ^ rijcr (2) 

{ij),cT i i,cr 

describes the correlated layer. {i,j) stands for neighbor¬ 
ing i and j sites, c| ^ creates an electron at the i-th site 

of the correlated layer with spin cr =t, i and njg. = 
denote the corresponding occupation-number operators. 
The leads are described by the Hamiltonian 

= —ta 4iia'^aja + 4iia^aia ' (^) 

Here creates an electron at i-th site of the lead a. An 
applied bias voltage $ shifts the energies £« and chemical 
potentials fXa of the leads in opposite directions by the 
amount $/2. Finally, 

"Hcoup = - Va (cLcaj> + h.C.^ (4) 

describes the hybridization between the correlated layer 
and leads. The hopping Va takes place between neigh¬ 
boring sites of the lead and the correlated layer. 

Previously similar models with many correlated layers 
were also investigated in Refs. |3Z1 [3S] [BH and . In 
Refs. IB71 and IBBl steady state behavior, while in Refs. [B¥l 
and[B^full time evolution were investigated. For this pur¬ 
pose the authors used DMFT IRefs. 15711551 and lB^ and 
time-dependent Gutzwiller approximation (Ref. I65II . In 
the Refs. [37l and [38] the impurity problem is treated by an 
equation-of-motion approach with a suitable decoupling 
scheme for the higher order Green’s functions, while in 
Ref. [BD the non-crossing approximation is invoked. On 
the other hand, our treatment of the impurity solver is 
controlled and can achieve extremely accurate result^^ 
with a moderate number of bath sites. 


A. Model 


B. Steady-state non-equilibrium Green’s functions 


To illustrate the approach we consider a minimalis- 
tic model for transport across a correlated interface (see 


We consider an initial situation in which at times r < 0 
the leads are disconnected from the correlated layer and 
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all three parts of the system (/, c, r) are in equilibrium 
with different values for the chemical potential /r; = £i, 
/ic = £c and = £r, respectively. 

Due to the fact that the system is transitionally invari¬ 
ant in the xy plane, it is more convenient to perform a 
Fourier transformation and express the Green’s functions 
in terms of the momentum k|| = {kx,ky). The retarded 
equilibrium Green’s function for the disconnected non¬ 
interacting central layer reads 




1 

LO + *0+ - £c- EcCkw) ’ 


(5) 


with i?c(k||) = —2tc{coskx + cosky). On the other hand, 
the Green’s functions for the edge layers of the left (a = 1) 
and the right (a = r) lead, when they are disconnected 
from the central layer can be expressed ad^^^H 




uj £(^ F/Q,(k||) 


-{Uj-£a- Fla(k||)P 


2tl 


( 6 ) 


with ifQ,(k||) = —2ta{coskx + cosky). The sign of the 
square-root for negative argument must be chosen such 
that the Green’s function has the correct l/w behavior 
for |a;| —>■ oo. To investigate the system out of the equilib¬ 
rium, we need to work with in the Keldysh Green’s func¬ 
tion formalism.SnilllilMIznixherefore, as a starting point, 
we need the corresponding non-interacting, disconnected 
Keldysh components. Since the disconnected systems are 
separately in equilibrium, we can obtain these from the 
retarded ones via the fluctuation dissipation theoreirP^ 


5f(w,k||) =2i(l-2/„(a;))Imgf(a;,k||). (7) 


Here, fa{ijj) is the Fermi distribution for chemical po¬ 
tential and temperature Tq,. For the non-interacting 
isolated central layer, the inverse Keldysh Green’s func¬ 
tion k||)]^ is infinitesimal and can be neglected 

in a steady state in which the layer is connected to the 
leads. In our notation, we use an underline to denote 
block matrices within the non-equilibrium Green’s func¬ 
tion (Keldysh) formalism: 


X = 


0 x^ j 


( 8 ) 


with X^ = (X^)I. At time r = 0 the leads get connected 
to the correlated layer. After a sufficiently long time a 
steady state is reached. The latter is expected to exist 
and to be unique unless the system has bound states. 
Our goal is to investigate its properties under the bias 
voltage $. 

Since the steady state is time-translation invariant, we 
can Fourier transform in time and express all Green’s 
functions in terms of a real frequency uj. The Green’s 
function for the correlated layer, when connected with 
the leads, can be expressed via Dyson’s equation 

G-^(a;,k||) = Go-'(a;,k||)-S(..,k||), (9) 
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FIG. 2. (Color online) Sketch of the auxiliary open quan¬ 
tum impurity problem consisting of the impurity site at po¬ 
sition i — 0 (red circle), Xj, = 6 bath sites (cyan circles) 
and two Markovian environments sink (blue) and reservoir 
(green). Parameters Eij and Fij are explained in the main 
text. 


where S(a;,k||) is the is self-energy of the correlated 
layer. The Green’s function of the non-interacting non¬ 
equilibrium system GQ(a;,k||) in turn can be expressed 
as 


Go^(w,k||) =£jji(a;,k||)- ^ vl gjuj^k^), (10) 

a—l,r 

where (/jj(a;,k||) is the Green’s function of the non¬ 
interacting decoupled layer, i.e. U = 0,Va = 0 and the 
components of the Green’s function of the isolated leads 
5 ^(a;,kj|) are given in Eq. (j^ and Eq. (0. Note that all 
quantities are underscored, i.e. they are Keldysh block 
matrices. 


C. Dynamical mean-field theory 

As usual, to obtain the self-energy E(a;, k||) is the dif¬ 
ficult step in the calculation of G(a;, k||) and of various 
steady state properties of the system. As there is no 
closed expression for it, one has to reso rt to some a p- 
proximation. Here, we employ DMFltlStSiElinSISIIlI] 
its nonequilibrium, time-independent version. In this ap¬ 
proach the self-energy is approximated by a local quan¬ 
tity S(a;, k||) = E(w) which can be determined by solving 
a (non-equilibrium) quantum impurity model with the 
same Hubbard interaction U and on-site energy £c cou¬ 
pled to a self-consistently determined bath. The latter is 
specified by its hybridization function obtained as 

A(a;) =£-!(..)-GfocH-Hc^), (11) 

where g~^{uj) is the non-interacting Green’s function of 
the disconnected impurity (i.e. of a single correlated site) 
and 


^G(w,k||). (12) 

BZ 

The self-consistent DMFT loop works similarly to the 
equilibrium case, except that in the pre sent cas e the 
Green’s functions are 2x2 block matricej ^^ * ^^*^^ ! One 
starts with an initial guess for the self-energy E(a;), then 
based on Eqs. calculates the bath hybridization 
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function A(w). We then evaluate the corresponding aux¬ 
iliary Green’s functions Ga^ux,o(^) Gaux(‘^)j ^^e 
non-interacting and in the interacting case, respectively. 
The solution of the impurity problem is, as usual, the 
bottleneck of DMFT. Our scheme consists, as outlined 
in detail in Sec. |IID[ in replacing the impurity problem 
with an auxiliary one, which is as close as possible to the 


one described by (111 but is exactly solvable by numer¬ 


ical methods. The self-consistent loop is then closed by 
determining the new value of the self-energy 


S(a;) — Gaux,o(‘^) “ Gaux(^) • 


(13) 


We repeat this procedure until convergence is reached, 
i.e. until Gaux(‘^) « Gioc(w)-‘^ 


D. Impurity solver: auxiliary master equation 
approach 

As already mentioned, the main obstacle of DMFT is 
the solution of the impurity problem. One widespread 
approach to approximate its solution for the equilibrium 
case is ED-DMFT, whereby one replaces the infinite bath 
with an auxiliary finite one. However, this approach can¬ 
not be used straightforwardly in a nonequilibrium steady- 
state case, as this cannot be described correctly with a 
finite number of sited^. One way to overcome this prob¬ 
lem, as some of us already suggested in Refs. [51] and [52] is 
to introduce, in addition to the finite number W of bath 
sites which are coupled to the impurity in the form of 
two chain segments, two Markovian environments, which 
can be seen as a particle sink and reservoir, respectively 
(see Fig. [^. This makes the impurity model effectively 
infinitely large, which is necessary in order to be able to 
reach a steady state. Our strategy, similar to the equi¬ 
librium ED-DMFT case, is to choose the parameters of 
the auxiliary model so as to provide an optimal fit to the 


bath hybridization function (111. 

The dynamics of this auxiliary impurity model, is de¬ 
scribed by the Lindblad quantum master equation, which 
controls the time (t) depen dence of the reduced density 
matrix p of the modeP^Ell 


^ /■ 
■^P = -Cp, 


where 


C = C 


H 


L 


D 


(14) 


(15) 


where d\^ creates a particle with spin a at the impurity 

[i = 0) or at a bath site (i = 1,..., W)- Jt-qo- = ^oo-'^Otr 
is the occupation number operator for particles at the 
impurity-site with spin a. Eqq = Sc, while all other Eij 
are parameters used to fit A(a;), whereby one can restrict 
to on-site and nearest neighbor (n.n.) terms only (see 
Fig.§. The non-unitary (dissipative) term 




CdP = 2 ^ ^ 

i,j —0 o 


.(1) 


(^iaP 


[d]^P - \{P^d^ad]a}^ 


(18) 


describes the coupling to a Markovian environment. The 
dissipation matrices k = 1,2 (with matrix elements 
r)^^) are Hermitian and positive semidefinit^^ and are 
again used as fit parameters. In order to fix the large-w 

behavior of A, all F,)^^ with at least one index on the 

- . 

impurity must vanish. On the other hand, in contrast to 
E, are not restricted to n.n. terms. This is of great 
advantage for the fit, as discussed below in Sec. [Ill A[ 

To carry out the self-consistent DMFT loop we need 
to evaluate both the non-interacting and the interacting 
Green’s functions of the auxiliary model. First, the non¬ 
interacting calculation (U = 0), which is fast in compar¬ 
ison to the interacting one, produces the bath hybridiza¬ 
tion function A„„y(a;) of the auxiliary impurity model. 


which is fitted to (111 in order to obtain the optimal pa¬ 
rameters E and These are used in the interacting 

model in order to determine the self-energy E(w), which 
is then inserted in 

A convenient way to solve the auxiliary problem is to 
rewrite Eq. (14|, expressed W superoperators, into a 
standard operator probleirp3HZ2l por this purpose one 
enlarges the original Fock space, spanned by the oper¬ 
ators (dj^/dj^), by doubling the number of levels via 
so-called tilde operators d^ „/d\^. In addition one intro¬ 
duces a so-called left vacuum 


\I) = ®\S). 


(19) 


where \S) are many body states of the original Fock 
space, \S) the corresponding ones of the tilde spac^^ 
and Ns the number of particles in S. In this formalism 
the reduced density operator is mapped onto the state 
vector 


is a Lindblad superoperator, which consist of two terms: 
a unitary contribution Cr and a dissipative one Ld. 
The unitary contribution 

ChP = -*[Haux, p] , (16) 

is generated by the auxiliary Hamiltonian 

Nb 

Tdau-K = ddijdj^dj^ + U Uq^Hq^ , (17) 


\p{r))=p\I), (20) 

and the Lindblad equation is mapped onto a Schrodinger- 
type equatiojJ^ 

±\p(r)) = L\p{r)), (21) 

where 

L = Lo + Li (22) 
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is an ordinary operator in the augmented space. Its non¬ 
interacting part Lq reads 

zLo = 5] (dt hd, - Tr(E + *A)) (23) 

cr 

where Tr denotes the matrix trace, and 

= (^ 0 .(T> ■ ■ • > ^ 0 , 0 -) • ■ • > ■ (24) 

is a vector of creation/annihilation operators and the ma¬ 
trix h is given by 


TVs > 3 due to the fact that the Hilbert space is expo¬ 
nentially large. Particle conservation translates here into 
conservation of — N^. To calculate the steady state 
I Poo) we use an Arnold! time evolutiorP^, while for the 
calculation of Green’s functions we employ the two-sided 
Lanczos algorithm.!^ 

To obtain the retarded and the Keldysh Green’s func¬ 
tions we use the following relations (for details see 
Ref. EU : 

Qi? ^ Q>+ ^ Q<-St 

= G>+ -h G<+ - H.c.. (34) 


h = 


E -h 2r(2) 
_2r(i) E-ift 


with 


A = r(2) -H , n = . 

Its interacting part has the form 

iLi = Uno^noi - Ufio^hoi, 


The expressions for the greater and lesser Green’s func- 

(25) 

tions ar^^ 




(26) 

and 

(27) 

g; 




_/(+!) 


(j|4|i?t'))(Lt'V,Jpoo) 


-I 


(-1) 


, (35) 


, (36) 


with fioo- := dlfadoa- evaluate Green’s functions, one 
needs to calculate expectation values of the form 

Gba =-R^u{B{T 2 )A{Ti)f>u{Ti)) , (28) 

where pu{ti) is the density operator of the “universe” U 
composed of the “system” (the chain in Fig. and the 
Markovian environment and tr = tr 0 tr ^ is the trace 
over the “universe”, which is the tensor product of the 
trace over the “system” (tr ) and the trace over the en¬ 
vironment (tr£)). After straightforward calculations we 
obtain (more details see Ref. for the non-interacting 
retarded Green’s function 

Gfux.o = (^-E + *A)-' (29) 

while its Keldysh part reads 


with the right (ji?),^^^)) and left ((l))^^^|) eigenstates and 
eigenvalues of the operator L (Eq. @) , in the 

sectors N„ — = ± 1 . 


Once self-consistency in the DMFT loop (cf. Sec. IIGI 
is achieved, one can calculate desired physical quantities, 
e.g. the steady-state current . For th is purpose we use 
the Meir-Wingreen expressiorP2l^2llI] jtg symmetrized 
form, where summation over spin is implicitly assumed. 


BZ 

+ (7z(w,k||) - >(a;,k||))(G^(a;,k||) - G"^(w,k||))] , (37) 


here 70 ,( 0 ;,k||) = - 2 u^Im 5 o(w,k||) and 70 ( 0 ;,k||) = 
f{uj - /ro) 7 o(a;,k||). 


G 


K 

aux,0 


= 2iG 


R 

aux,0 


QG 


A 

aux,0 ’ 


(30) 


III. RESULTS 


Therefore, we obtain the following expressions for the 
retarded auxiliary hybridization function 


Afux = w - Ec - pTT 


,oJoo 


and its Keldysh part 


^aux 


= 


[Gaux,o]oO 


G^x.olooP 


(31) 


(32) 


Here Ago denotes the 00 element (i.e. the one on the 
impurity) of the matrix X. 

To calculate the impurity Green’s function for the in¬ 
teracting system we use Krylov-space based exact di- 
agonalization. A full diagonalization is prohibitive for 


Here we present results for the steady state properties 
of the system displayed in Fig. consisting of a corre¬ 
lated layer, with Hubbard interaction U and on-site en¬ 
ergy Ec = —U/2, coupled to two metallic leads. We re¬ 
strict to the particle-hole symmetric case. The hopping 
inside the correlated layer tc is taken equal to 1 unless 
stated otherwise, while the hopping inside the leads is 
U = tr = 2.5. The hybridizations between leads and 
the correlated layer are vi = Vr = 0.5 and 2vi is used as 
unit of energy. The applied bias voltage 4> enters the val¬ 
ues of the onsite energies and the chemical potentials as 
Ei/r = pi/r = ±$/2. All results presented below are cal¬ 
culated for zero temperature in the leads (T; = Tr = 0 ), 
see Eq. (Q. Similar models have been studied, e.g. in 
Refs. I37n381 1551 EH and ESI 
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FIG. 3. (Color online) Retarded and Keldysh Green’s 
functions for the correlated layer, for Hubbard interaction 
U = 2 and bias voltage $ = 5 (upper panel) and for U = 12 
and $ = 10 (lower panel). Other parameters are: 1^ = 1, 
ti = tr = 2.5, vi = Vr = 0.5 {2vi is our unit of energy). Solid, 
dashed and dotted lines are obtained with A^i, = 6, A^;, = 4 
and Ni, = 2 correspondingly. 


A. Convergence with respect to the number of 
auxiliary bath sites Ni, 

First we investigate how the number of bath sites Nb of 
the auxiliary impurity problem influences the results. We 
compare calculations for the Green’s functions (Fig. 
and for the current (Fig. |^, obtained with Nb = 2,4,6. 
We And that the retarded component is well converged 
already for Nb = 4 even for U = 12. For the Keldysh 
Green’s functions the convergence in terms of Nb is rea¬ 
sonable, but not as fast as for the retarded Green’s func¬ 
tion. Gorrespondingly, it is not surprising that also the 
current voltage characteristics exhibit a fairly good con¬ 
vergence (See Fig. |^. On the whole, the convergence 
for weak interaction (U = 2) is faster than for strong 
interaction {U = 12). 

These results indicate that Nb = 4 bath sites already 
produce reasonable results away from the Kondo regime. 
Therefore, in view of the exponential increase of the nu¬ 
merical effort with Nb, we mainly restrict the following 
discussion to Nb = 4. Only to discuss the low energy 
Kondo physics we will present results with Nb = 6. The 
reason for the rapid convergence in Nb is due to the 
fact that the number of Lindblad parameters increases 
quadratically with Nb, in contrast to the energy and hop¬ 



FIG. 4. (Golor online) Gurrent J vs bias voltage for dif¬ 
ferent U. The solid red curve represents the U = 0 result. 
The three curves peaked at about il? = 4 are for U = 2 and 
the remaining three curves show the 17 = 12 results. Solid, 
dashed and dotted lines correspond to Nb = 6, Nb = 4 and 
Nb = 2, respectively. Other parameters are as in Fig. 



FIG. 5. (Golor online) Gurrent J vs interaction U for two 
values of 4>. Galculations are performed for Nb = 4. Other 
parameters are as in Fig. 

ping parameters E. It is, thus, important, to consider 
also long-ranged T terms (cf. Ref. 15^ 


B. Steady-state current 

In this section we discuss the steady-state current in 
detail. The results for the current as a function of bias 
voltage are presented in Fig for different values of 
the Hubbard interaction U. In the non-interacting case 
{U = 0) particles pass the interface without scattering 
and therefore the momentum k|| is conserved. Gorre¬ 
spondingly, the problem becomes one-dimensional and 
the current vanishes for bias voltages larger than the one¬ 
dimensional bandwidth, i.e. for HA, = 4ti/r = 10, which 
is corroborated by Fig. For nonzero interaction dif¬ 
ferent k|| are mixed due to scattering and thus all states 
of the leads are possible Anal states. Subsequently, the 
current vanishes for bias voltages larger than the three- 
dimensional bandwidth, i.e. W = 3 ■ Wz = 30. In equi- 
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librium, an isolated two-dimensional Hubbard layer is in 
the metallic phase for weak interaction. As can be seen 
from Fig. in this case {U = 2) the current displays, as 
expected, a metallic behavior, i.e. a linear increase of the 
current for small voltages. The overall shape is similar to 
the U = 0 case, however, with a longer tail at large $ due 
to the scattering mechanism discussed above. For strong 
interaction {U = 12) an isolated two-dimensional Hub¬ 
bard layer is a Mott insulator, but in our model there 
is no insulating phase due to the hybridization to the 
non-interacting leads. Therefore, strictly speaking the 
current is always linear in d) for $ —>■ 0. Nevertheless, 
due to the vicinity of the Mott insulator the current is 
strongly suppressed.^ A similar behavior also was ob¬ 
served in Refs. 1571 and 155] . On the other hand, for higher 
bias voltages ($ > 12) the picture is reversed and the 
current is more suppressed for 17 = 2 than for U = 12. 

We investigate this issue in detail and plot the current 
as a function of Hubbard interaction U for low ($ = 3.5) 
and high ($ = 15) bias voltage in Fig. For the low 
voltage case, we find a monotonic Gaussian decrease. 
The origin of the current reduction with increasing U are 
back-scattering processes that reduce the transmission 
coefficient. For high bias voltage, in the region where 
the current is zero for U — 0, the current first increases 
with increasing interaction, reaches its maximum at ap¬ 
proximately 17 ~ 6 and then decreases again.l^ Qual¬ 
itatively this can be explained by the fact that there 
are two competing effects as a function of U. On the 
one hand, with increasing U the transport increases due 
to scattering to different k|| as discussed above, which 
enhances the current, but on the other hand, large U 
means increased backscattering which suppresses trans¬ 
port across the correlated layer. For high bias voltages 
and weak interactions the first effect dominates due to 
the finite bandwidth of the leads. 



FIG. 7. (Color online) Single particle spectral function for 
Ni, = 6, U = 2, tc = 0.1 and for different values of $. Other 
parameters are as in Fig. 



FIG. 8. (Color online) Single particle spectral function for 
Nb = 6, U = 12, tc = 1 and for different values of $. Other 
parameters are as in Fig. 


C. Non-equilibrium spectral function 

To gain further insight into the properties of the steady 
state we also investigate the non-equilibrium spectral 
function, which can be calculated from the Green’s func¬ 
tion via A{uj) = The results are shown in 


















(a) (b) (c) 

FIG. 9. (Color online) Single particle spectral function for Nb = 4, different values of 17, 4> = 0 (a), 4> = 3.5 (b) and >1> = 15 
(c). Other parameters are as in Fig. 


Fig.i for 17 = 2 and U = 12. For weak interaction 
(17 = 2) the spectral fnnction A{uj) displays for all bias 
voltages a peak at w = 0 and hardly visible Hnbbard 
satellites at the approximate position ui = zLUl2. The 
spectral function for U = 2 (Fig. 6(a)) depends only 
very weakly on bias voltage. This is in contrast to the 
spectral density of the one-dimensional SIAM for which, 
it is found that, the Kondo peak splits up as a function 
of voltage and two resonances are observed at the cor re- 
sponding chemical potentials of the two lead j'^^F5 | 62 | 86| (9l] 
The difference between the two situations is that in the 
case of the single impurity model the resonance and the 
Hubbard subbands are clearly separated. In contrast, for 
the correlated layer and for the set of parameters consid¬ 
ered here, when the isolated correlated layer is metallic, 
they overlap due to the broadening induced by the hop¬ 
ping tc within the correlated layer. And indeed if we 
artificially reduce tc to 0.1 (keeping all other parameters 
fixed) we observe, for $ = 0 a resonance, which is clearly 
separated from the Hubbard subbands (cf. Fig. [^. In 
this case, the isolated layer would be insulating and the 
broadening of the resonance is not any more due to tc, 
but to an effective energy scale Tk, which can be seen as a 
Kondo temperature. This originates from a combination 
of coherent scattering from the leads into the insulating 
layer, as well as from a self-consistent DMFT process 
as discussed in Ref. [13 In addition to the broadening 
mentioned above there is also a spurious broadening due 
to the limited accuracy of our calculation. Nevertheless, 
our resolution is sufficient in order to observe a splitting 
of the resonance into two peaks at ^i/j. = ±.Vb/2. as a 
function of voltage as in the single impurity case. 

Now we turn to strong interactions {U = 12), for which 
the results are depicted in Figs. |6(b)| and In equilib¬ 
rium, i.e. for $ = 0, the hybridization with the leads 
produces a weak Kondo resonance at the Fermi Energy 
(w = A nonzero bias voltage splits the resonance 

into two peaks at = ±$/2 (see Fig.l^. For $ > 3 the 
peaks merge into the Hubbard subbaiid^s and the spec¬ 
tral function A{uj) consists of two Hubbard subbands at 
the approximate position oj = ±U/2, while A(a; = 0) is 
strongly suppressed. For these larger values of $ the ef¬ 


fect of the bias voltage is small: it modifies only slightly 
the position and height of the Hubbard subbands. No¬ 
tice that in order to resolve the Kondo peak and its split¬ 
ting at low bias we need to use an auxiliary system with 
Nb = 6. While this allows to resolve the peaks, the lim¬ 
ited accuracy makes them broader than they should be 
and therefore the Friedel sum rule is not satisfied even for 
equilibrium. The reason for this is that a spurious broad¬ 
ening originating from the limitation of our approach re¬ 
duces the height of the peak at the Fermi level.l^To fulfill 
this one would have to use more bath sites, and conse¬ 
quently adopt a matrix-product state based solution of 
the auxiliary system, as we did in Ref. [53] for the Ander¬ 
son impurity model. This, however, in combination with 
the DMFT self-consistency, would increase considerably 
the required computation time. 

Next we study the dependence on the Hubbard inter¬ 
action U in more detail. Results for equilibrium ($ = 0), 
low ($ = 3.5) and high ($ = 15) bias voltages are pre¬ 
sented in Fig. As one can see, in the equilibrium case 
for large U only small excitations are visible at a; = 0. 
Here the Friedel sum rule would require the w = 0 peak 
to display a [7-independent height.!^ However, as dis¬ 
cussed above for large U our approach cannot resolve Tk 
and the peak becomes strongly suppressed, and for in¬ 
termediate U the sum rule is not satisfied, as for Fig. 
Still, Fig.|^a) displays a crossover from a regime in which 
the local Fermi liquid peak is already present in the iso¬ 
lated correlated layer (for U below the 2D Mott tran¬ 
sition), into a Kondo-Fermi liquid regime in which the 
peak is produced by coherent spin-flip processes across 
the Mott insulator, originating from the Fermi levels of 
the leads. Outside of the Kondo regime, the behavior 
is qualitatively similar in the three cases. In the non¬ 
equilibrium cases, upon increasing the interaction the 
height of the spectral function at w = 0 decreases and 
for [7 > 6 the spectral function displays a local minimum 
at w = 0 instead of a maximum, which becomes vanish¬ 
ingly small with increasing U. Comparison of Figs. |9(b)| 
and |9(c)| shows that for higher bias voltages this res¬ 
onance disappears at smaller U. Also, for higher bias 
voltages (<!> > 10) the non-interacting spectral function 
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has a sharper peak at w = 0. This is due to the fact 
that for high bias voltage the leads’ density of states do 
not overlap any more and correspondingly states close to 
the Fermi level do not dissipate any more into the leads, 
therefore the density of states close to the Fermi level is 
just the two dimensional density of states, which features 
a logarithmic divergence. Another effect, clearly visible 
in Figs. 9(a) 9(b) and 9(c)[ is the linear shift of the po¬ 
sition of the Hubbard subbands with increasing U. The 
peak position is given by uj± ~ ±17/2. 


IV. CONCLUSIONS 


We have presented an improved application of a DMFT 
technique for non-equilibrium situations that allows to 
study directly steady state properties of strongly corre¬ 
lated devices. Like in equilibrium DMFT, the only ap¬ 
proximation is the locality of the self-energy, while the 
accuracy of the non-equilibrium impurity problem is con¬ 
trolled by the number of bath sites W which are attached 
to Lindblad environments. We find that the accuracy in¬ 
creases exponentially with both in and out of equi¬ 
librium. The approach is benchmarked for a strongly 
correlated layer coupled to two metallic leads. While the 
results in Ref. EB for this model were obtained by full 
diagonalisation of the auxiliary impurity problem and 
were, thus, restricted to Vf, = 2, here we invoked the 
non-hermitian Krylov-space method, which allows us to 
use larger values for N^. 

With the Krylov-space solver we were able to go up to 
Nb = 6. For more bath sides (up to Nf, ^ 14) the MPS 
solvei)^ could be used, but the Krylov-space solver has 
the advantage to be quicker and more flexible, which is 
important for the DMFT iteration. For the single layer 
device studied here we found that Aj, = 4 already yields 


very reliable results in most parameter cases. Only the 
Kondo regime requires larger values for W, but with 
fVb = 6 at least semi-quantitative results can be achieved. 

We have investigated the current-voltage characteris¬ 
tics across a correlated layer. At low bias voltages, we 
have observed a linear behavior for weak interactions, 
while the current was exponentially suppressed for strong 
interactions.!^ On the other hand, for higher bias volt¬ 
ages we have observed a reversed picture, whereby the 
current is larger in the strongly interacting case. In ad¬ 
dition we have investigated the current J as a function of 
the local Hubbard interaction for low as well as for high 
bias voltages. For lower bias voltages we found that the 
J decreases monotonically with U, while for higher bias 
voltages, the current first increases, reaches its maximum 
and then decreases again. The origin of this behavior can 
be explained by different scattering processes. 

In addition to the current we have also investigated 
the steady state spectral function. Our results show that 
for the set of parameters considered in this manuscript, 
the spectral function is only weakly dependent on the 
bias voltage in contrast to the single impurity problem. 
This is due to the fact that the splitting of the Kondo 
resonance as a function of $ is strongly broadened due 
to the hopping within the correlated layer. As to be 
expected, the Hubbard satellites depend almost linearly 
on U, like in the equilibrium case. 
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